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ON  THE  CONVERGENCE  OF  THE  CONJUGATE  GRADIENT  METHOD 


FOR  SINGULAR  CAPACITANCE  MATRIX  EQUATIONS* 
A.  S.  L.  Shieh 


§ 1 . Introduction 

Over  the  past  decade,  very  fast  direct  methods  have  been  developed  to  solve  Poisson's 

equation  on  certain  simple  regions  with  Dirichlet,  Neumann  or  periodic  boundary  conditions. 

See  c.g.  [2],  [7],  [lO],  [l7],  ll9],  and  [ 29 ) • The  capacitance  matrix  methods  are 

developed  recently  to  solve  Poisson's  equation  on  arbitrary  bounded  regions  with  smooth 

boundaries  by  imbedding  the  discrete  problem  into  a region  where  these  fast  direct  methods 

are  applicable.  See  section  7 of  this  work  for  a brief  survey  of  previous  work  on  capacitance 

matrix  methods.  In  this  work  it  is  shown  mathematically  that  by  making  the  correct 

Ansatz  guided  by  classical  potential  theory,  the  convergence  of  the  conjugate  gradient 

method  for  solving  the  capacitance  matrix  equations  is  essentially  independent  of  the  mesh 

2 2 

size.  The  total  operation  counts  of  the  algorithm  do  not  exceed  constant  n (log  n)  where 
h = l/n  Is  the  mesh  size.  Only  numerical  schemes  of  first  order  accuracy  for  the  interior 
Neumann  problem  of  the  Poisson  equation  on  bounded  two  dimensional  regions  with  smooth 
boundaries  arc  considered  here.  See  [ 28)  for  a similar  treatment  of  the  Dirichlet  problem. 
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We  give  only  a very  brief  review  of  a few  results  of  classical  potential  theory.  For  a 


i 

i 

I 


detailed  exposition  see  e.g.  [8],  [l2),  [22]  and  [25).  Wo  define  the  potential  1r 
resulting  from  a charge  distribution  p on  a smooth  boundary  curve  dU  by 

r(x)  = (1/tr)  f p(4)  log  r ds(|)  . 

an 

Here  x = (x^,x^),  i = (ij,  and  r^  = (x^  - 4^)^  + (x^  - The  Green's  function 

(l/2n)  log  r which  we  shall  denote  by  G satisfies 

A(1/2tt)  log  r = 6{r)  , 

where  6(r)  is  the  delta  function.  For  the  interior  Neumann  problem,  we  make  the  Ansatz 

(2.1)  u(x)  = (1/2TT)  //f(4)G*d|  +?-(x) 

n 


for  the  solution  of 


(2.2) 


Au  = f,  X t J2 


au/dv  = g,  X < an  . 

Hero  V denotes  the  outer  normal  to  the  boundary  curve  an.  The  first  term  on  the  right 
hand  side  of  (2.1)  is  a space  potential  term  and  will  be  denoted  by  u„.  The  boundary 

O 

condition  is  satisfied  by  choosing  p such  that 


(2.3)  p - (1/ir)  J paG  /av  ds  = g + (a/av)u  I 

an  ^ ^'an 

This  equation  can  be  written  as 


(2.4)  (1-K)p  = 5, 

where  K is  a compact  operator  defined  by  the  integral  above.  The  equation  Is  a Fredholm 

integral  equation  of  the  second  kind  and  thus  a well  posed  problem.  It  has  a simple  zero 

~ * 

eigenvalue  and  is  solvable  if  g has  a zero  mean  value.  We  remark  that  G in  equations 

(2.1)  and  (2.  3)  can  be  replaced  by  the  Green's  function  on  a rectangle  with  zero  Dirichlet 
boundary  conditions  or  any  other  Green's  function  of  the  Laplacian. 
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§ 3.  The  Ccipnritance  matrix  method  for  the  Neumann  problem 

In  this  section  we  develop  a similar  formal  potential  theory  for  the  discrete  problems 
arising  from  the  original  Neumann  problem  (2.2).  See  also  sections  3 and  5 of  [26]  for  a 
similar  discussion.  We  shall  assume  that  uniform  mesh  sizes  in  both  coordinate  directions 
are  used. 

* 

We  replace  the  Laplace  operator  by  the  five-point  formula.  The  Green's  function  G 

used  in  section  2 will  then  be  replaced  by  the  discrete  Green's  function  on  a rectangular 

region  S with  Dirichlet  boundary  conditions.  We  denote  this  Green's  function  by  B ^ 

2 

where  B is  the  matrix  representing  the  discrete  Laplacian  h A^,  employing  undivided 

differences,  on  S and  zero  boundary  values  on  the  grid  points  of  9S. 

We  imbed  U in  S as  follows.  The  set  of  mesh  points  is  decomposed  into  three 

disjoint  sets  n.  , dU.  and  (CO)  . The  set  gO.  contains  all  the  Irregular  mesh  points 
h h h h 

in  n,  i.e.  mesh  points  that  do  not  have  all  four  neighbors  within  the  open  set  O. 

is  the  sot  of  regular  mesh  points  inside  U and  (CO,  ) contains  the  remaining,  the 

h 

exterior  mesh  points.  We  further  require  that  O is  bounded  away  from  aS  uniformly  in 

h.  We  then  set  up  tiro  matrix  equation 

(3. 1)  Au  = V 

that  we  are  solving  as  follows.  We  require  that  B and  A differ  only  on  the  rows  that 

corresponds  to  the  irregular  mesh  points.  On  these  rows  we  combine  the  discrete 

Laplacian  with  difference  approxinjatlons  to  the  normal  derivative.  We  must,  however,  be 

sure  that  the  solution  on  n.  U an.  is  Independent  of  the  solution  or  data  on  (Cn),  . 

h h n 

This  is  achieved  by  eliminating  from  the  discrete  Laplacian,  centered  at  an  Irregular  mesh 
point,  the  values  of  the  solution  at  its  exterior  neighbors.  We  write 

A = B - UV^ . 

/ 
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The  matrices  U and  V have  m columns  where  m is  the  number  of  points  in  . 

h 

U represents  an  extension  operator  that  retains  the  values  of  mesh  function  on  dn.  and 

h 

T 

makes  the  remaining  values  equal  to  zero.  The  rows  of  V are  simply  the  differences 
between  the  corresponding  rows  of  B and  A.  After  a suitable  permutation,  the  matrix 
A is  reducible, 


(3.2) 


The  submatrix  A^  ^ is  the  matrix  for  the  linear  system  of  equations  of  the  original  discrete 

problem  arising  from  discretization  of  the  original  problem  (2.2).  It  is  easily  seen  that 

the  restriction  to  n.  U an.  of  any  solution  of  Au  = v must  be  a solution  to  the  original 
h h 

T T 

discrete  problem.  V will  be  chosen  so  that  the  row  sums  of  A^  ^ and  V vanish  and 

Aj  j has  a simple  zero  eigenvalue.  The  matrix  A^^  is  nonsingular  since  it  represents 

a finite  difference  approximation  to  a Dirichlet  probiem  on  CQ.  It  is  then  easily  verified 

that  the  matrix  A also  have  a simple  zero  eigenvalue. 

We  now  describe  our  method  for  solving  the  system  equations  (3.1).  It  is  solvable 

If  and  only  if  the  right-hand  side  v is  orthogonal  to  the  left  eigenvector  of  A which 

corresponds  to  the  zero  eigenvalue.  It  is  shown  in  section  5 of  [ 26]  that  the  right  hand 

side  V is  always  consistent  regardless  of  its  values  on  (Cfl^^)  if  the  data  is  already 

consistent  on  il,  U ail,  . 

h h 

Guided  by  the  classical  potential  theory,  we  make  the-  Ansatz 


(3.  3)  u = b"'v  + b"'  UDp  . 

Here  p is  a m-vector  to  be  determined.  D is  a nonsingular  diagonal  matrix  containing 


certain  scaling  factors  to  be  specified  later.  Computing  the  residual  vector  we  obtain 


Because  of  the  factor  U the  residuals  are  zero  for  all  x t They  must  also 

T 

vanish  on  We  therefore  multiply  equation  (3.  4)  by  U and  obtain 

(3.5)  (D  - V^b"'uD)p  = V^B~'v  . 

T 

Here  we  have  used  the  relation  U U = I , the  m x m identity  matrix.  We  shall  refer 

m 

to  equation  (3.  5)  as  the  capacitance  matrix  equation  and  the  matrix  on  the  left  hand  side 

of  (3.  5)  as  the  capacitance  matrix.  It  is  shown  in  section  5 of  [ 26]  that  the  capacitance 

matrix  which  we  shall  denote  by  C has  a simple  zero  eicjenvalue  and  that  tha  right  hand 
T - 1 

side  V B v of  (3.  5)  is  consistent  if  v is  consistent  for  the  original  problem  Au  = v. 

T - 1 

For  the  special  case  when  v = DU  v,  wo  can  simply  make  the  Ansatz  u = B UDp.  The 

capacitance  matrix  equation  now  becomes 

(3.6)  Cp  = u’^v  . 

T T TTT-l  TT-I 

Let  <()  satisfy  <)>  C = 0.  Then  4>  = <*>  V B U.  Therefore  <{)  V B v = 0 implies 

T T T 

()>  If  V = 0.  Hence  the  right  hand  side  U v of  equation  (3.6)  is  again  consistent  if  v is 
consistent  for  Au  = v. 

We  now  describe  our  choices  of  difference  equations  at  the  irregular  mesh  points. 

* 

Let  P t Let  P'  be  its  closest  point  on  aU.  Let  W,  E,  N,  S and  NE 

be  its  western,  eastern,  northern,  southern  and  northeastern  neighbors  on  the  mesh. 

Assume  that  the  local  orientation  of  the  boundary  is  such  that  W is  always  in 

while  N is  cither  in  gW.  or  (Ctl).  depending  on  whether  P has  one  or  two  neighbors 

n h 

in  (Cn)^.  Let  o <'n'/4  be  me  angle  that  the  normal  through  P makes  with  the  closest 
coordinate  axis.  We  approximate  the  Neumann  boundary  conditions  by  the  following  first 
order  scheme. 

(3.7)  u(W)  - (1  - tan  o)u(P)  - (tan  o)u(S)  = g(P  )h  cos  a 

(3.'8)  u(N)  - (1  - tan  n-)u(NE)  - (tan  a)u(E)  = g(P  )h  cos  a . 
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In  our  first  scheme,  we  combine  the  equations  (3.7)  and  (3.8)  with  the  discrete 


Laplacian  to  form  the  following  equation  regardless  of  whether  N is  in  (Cn).  or  not 

h 

(3.9)  cos  «[(3  + tan  o)u(P)  - (1  + tan  a)u(S)  - (1  - tan  ar)u(NL)  - (1  + tan  o)u(E)] 

= cos  o[h^f(I’)  + 2g(P  )h  cos  or]  . 

We  shall  refer  to  this  scheme  as  scheme  l.N.a. 

The  second  scheme  is  as  follows.  If  P € has  two  neighbors  in  (Cil)^^,  we 

obtain  the  equation  (3.9)  as  in  scheme  l.N.a.  If  only  W is  in  (Cn)^,  we  only  use 
equation  (3.7)  to  combine  with  the  discrete  Laplacian.  We  then  multiply  both  sides  of 
the  combined  equation  by  2 cos  a instead  of  cos  a to  obtain  the  following  equation 

(3.10)  2 cos  a[(3  + tan  o)u(P)  - u(N)  - u(E)  - (1  + tan  a)u(S)] 

* 

= 2 cos  a[f(P)  + g(P  )h  cos  a)  . 

We  shall  refer  to  this  scheme  as  scheme  I.N.b.  Both  schemes  l.N.a  and  I.N.b  give  rise 
to  matrices  that  are  positive  semidefinite  with  null  space  of  dimension  one  that 

consists  only  of  constant  functions.  It  is  easily  verified  (see  e.g.  [ 3))  that  the  solutions 
of  the  discrete  jiroblems  are  0(h  log  h)  approximations  to  the  exact  solutions. 

The  matrix  D on  the  left  hand  side  of  equation  (3.  5)  contains  the  scaling  factors 
dp  = sec  a.  Here  dp  is  the  diagonal  element  of  D on  the  row  corresponding  to  the  point 
P ‘ scaling  factors  cos  a and  2 cos  a in  equation  (3.9)  and  (3.10)  respectively 

and  the  diagonal  elements  dp  of  D are  chosen  so  that  the  off  diagonal  part  of  C may 

be  a formal  approximation  to  the  compact  integral  operator  K defined  by  equations 

* 

(2.3)-(2.4)  with  G replaced  by  the  Green's  function  on  a rectangle  with  Dirichlet  boundary 
conditions.  Because  of  the  irregular  patterns  of  points  in  near  diagonal  part, 

the  remaining  part  of  C,  will  not  in  general  be  a formal  approximation  to  the  identity 
operator.  It  will  be  shown  in  section  S that  this  near  diagonal  part  is  uniformly  well 
conditioned  in  the  spectral  norm  and  that  the  singular  values  of  C are  distributed  like  that 


t 

-6- 


of  the  sum  of  a positive  definite  symmetric  operator  and  a compact  symmetric  operator.  It 
is  shown  in  [14]  that  the  convergence  of  conjugate  gradient  method  for  solving  operator 
equations  with  such  operators  depend  asymptotically  only  on  the  spectral  condition  number 
of  the  positive  definite  symmetric  parts  of  the  operators.  The  method  of  proof  in  [14]  does 
not  apply  in  our  case.  We  shall  show,  however,  in  section  6 with  a different  approach 
that  the  corresponding  rates  of  convergence  in  our  cases  depend  also  asymptotically  on 
the  spectral  condition  numbers  of  B . 
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5 4 . Computationol  procotlures  and  oBeratlon  counts 


We  shall  solve  the  matrix  equation 

T T T - 1 

(4.1)  C Cp  = C V B V 

by  the  conjugate  gradient  method.  This  is  equivalent  to  solving  the  least  square  problem 

for  the  capacitance  matrix  equation.  See  section  6 for  further  details. 

T - 1 

In  principle  we  can  set  up  the  matrix  C by  computing  V B UD.  This  takes  at 

2 

least  m fast  Poisson  solvers  and  m storage  requirement.  It  is  therefore  much  better 

to  use  the  following  algorithm.  For  any  vector  v we  compute  Cv  as  follov/s.  Generate 

-1  T - 1 

the  mesh  function  UDv,  use  the  fast  solver  to  obtain  B UDv  and  compute  V B UDv 

T 

at  an  expense  on  the  order  of  m operations.  The  vector  C Cv  can  be  obtained  in  this 
fashion  at  a cost  of  essentially  two  fast  solvers.  It  is  easily  seen  from  (6.1)  and  (6.  2)  of 
section  6 that  each  Iteration  of  the  conjugate  gradient  method  will  therefore  cost  about  two 
fast  solvers.  The  theory  presented  in  sections  5 and  6 docs  not  preclude  the  possibility 
that  the  number  of  iterations  to  achieve  a given  accuracy  grows  like  log  m as  we  refine 
the  mesh  size.  We  have,  however,  consistently  found  in  our  experiments  that  the  number 
of  iterations  stays  constant  as  m increases  and  that  we  can  achieve  an  accuracy  of  between 
two  and  three  correct  decimal  digits  for  only  four  iterations.  The  operation  counts  for  many 
discrete  problems  are  therefore  ten  times  that  of  a fast  Poisson  solver  on  a rectangle  and 
the  storage  requirements  are  of  the  order  n^  whore  n = 1/h.  We  have  used  the  generalized 
marching  algorithm  described  in  [ 2 ] for  our  fast  solver  on  the  rectangle.  The  operation 
counts  of  this  fast  solver  is  approximately  3n^log2(~^j^)>  where  k is  the  size  of  each 
marching  step,  if  n = k2^  - 1 for  some  positive  integer  /.  The  marching  algorithm  is 
unstable  for  large  k.  Wo  have,  however,  found  that  k = 16  is  good  enough  for  our 
purpose.  The  operation  counts  for  many  problems  therefore  do  not  exceed  120  n^.  It  Is 
possible  that  if  fast  Fourier  transform  methods  are  used  to  compute  B *UDv  in  the 


computation  of  Cv,  the  algorithm  will  be  even  more  efficient  if  we  exploit  the  sparsity 
of  the  vector  UDv  in  the  Fourier  analysis  step.  One  big  advantage  of  the  capacitance 
matrix  method  is  that  it  can  be  speeded  up  by  the  replacement  of  a subroutine,  whenever  a 
faster  Poisson  solver  becomes  available.  Finally  we  mention  that  our  algorithm  can  be 
used  for  the  numerical  solution  of  the  Neumann  problem  for 

-Au  4 Cu  = f on  n,  C > 0 , 

although  the  theoretical  results  in  this  work  does  not  immediately  apply. 


•J 


§ 5.  The  distribution  of  singular  values  of  C 


We  shall  show  that  given  e > 0,  then  almost  all  the  singular  values  of  C lie  in 
the  interval  (dj  - c,  ^2  ^ where  dj  and  d^  are  positive  numbers  independent  of  h. 
This  is  accomplished  by  first  proving  that  is  uniformly  well  conditioned  in  the  spectral 
norm  and  that  almost  all  the  singular  values  of  lie  in  the  interval  [O,  e].  Our  main 

result  then  follows  by  a simple  application  of  a well  known  result  in  matrix  theory  which 
vjc  shall  state  below  as  lemma  5.8. 

IX^finition.  The  matrix  which  we  refer  to  as  the  near  diagonal  part  of  the  capacitance 

s 

matrix  is  defined  as  follows.  Each  entry  of  B^  that  corresponds  to  the  irregular  mesh 

h 

points  P and  O is  zero  if  d(P,  Q)  > s/^li ; otherwise  Bj^(P,  0)  hC(P,  O).  Here  d(  • , •) 
denotes  the  Euclidean  distance  function. 

Definition.  The  matrix  which  we  refer  to  as  the  off  diagonal  part  of  the  capacitance 

matrix  is  defined  to  bo  the  difference  between  M and  B,  , i.c.  K,  = C - B,  . 

h’  h h 

Theorem  5.1.  Given  e > 0,  there  exists  a positive  integer  N such  that  for  all  0 < h, 
all  except  N singular  values  of  lie  in  [0,e]. 

This  theorem  will  bo  a consequence  of  lemmas  5.1-5.  6 below.  First  we  need  some 
basic  results  from  the  theory  of  collectively  compact  operators.  Lot  K : X - X be  a compact 
operator  on  a complex  Banach  space  X. 

Definition.  A subset  S (-  X is  sequentially  compact  if  any  sequence  in  S contains  a 
convergent  subsequence  with  limit  in  X. 

Definition.  A family  of  operators  on  X is  collectively  compact  If  the  set 

{K^f  : Ilf  II  < 1,  f < X,  n = 1,  2, . . . ) is  sequentially  compact  in  X.  The  following  result 

is  an  immediate  consequence  of  a theorem  in  (1).  See  also  Chapter  4 of  [26]. 

Lemma  5. 1.  Let  (1^^)  3 family  of  collectively  compact  operators  on  a complex  Banach 

space  X with  converging  pointwlso  to  a compact  operator  K.  Given  e > 0,  let  pj, 
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with  algebraic  multiplicities  m.,  be  the  N(c)  eigenvalues  of  K with  absolute  values 

* * * 
greater  than  t.  Then  there  exists  a positive  integer  N and  an  t >0  such  that  for  all  n>N  , 

each  r neighborhood  of  contains  exactly  eigenvalues  of  while  all  the 

other  eigenvalues  of  K lie  in  an  e-neighborhood  of  zero. 

n 

We  now  construct  a family  of  operators  from  tis  follows.  Let  6,  ij) 

both  mappings  from  the  unit  interval  [0,1]  to  the  real  line,  be  a smooth  parametc  'zation 
of  Let  P.,  j = 1 m be  the  irregular  mesh  points . Let  P.  = (<>{t.),  i}j(t.))  be 

•k 

the  points  on  9f2  which  lie  on  the  normals  through  P^  with  d(P^,P^)  < h.  We  require  that 


0<tj<t^< 


< I < 1 . 
m 


Let  L denote  the  space  of  m tuples  v.-ith  the  sup  norm.  Let  C[0,1]  denote  the  Banach 
m 

space  of  continuous  functions  on  the  unit  interval  v/ith  the  sup  norm.  Let 


(5.1) 

where 

(5.2) 

(5.3) 

(5.4) 


Mt,t.)  ^ vV‘’i-r  ^ ' 

V = (t  - t^)/(t._j  - t.),  t._,  < t < t.,  i = 2,  . . .,m  : 

V « (I  - - 1),  t < I',,,  . P,  . P, 


P^.P, 


Define  bounded  linear  operators  P : C[  0 1 ] -*  L , K : L -♦  C|  0, 1 

m m rn  rn 

K : C[0,t]  - Cl0,l]  by 
m 


and 


(5.5) 

(5.6) 

(5.7) 


P^f=V,  V,  = 


m 

(K^v)(t)  r:  ^ k(t,t.)v(t.)  : 


)=1 


T 

We  then  construct  another  family  of  operators  from  above  procedure. 
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Let  K - K'  K . It  is  easily  verified  that  {K  ) is  also  the  family  of  operators  formed  from 
s m m s 

T 

{K.  K,  } by  the  same  procedure, 
h h 

Lemma  S.2.  Let  [X]  denote  the  Banach  space  of  bounded  linear  operators  on  a Banach 

space  X.  For  any  \ * 0, 

(M-K\)'*e[L  1 iff  (XI  - K r*elC[0,l]  ] . 
ii  n m s 

Proof.  For  X = 1,  the  lemma  is  proved  in  [ 1]  . Exactly  the  same  argument  applie.s  for  any  X ? 0. 

T 

Therefore,  the  nonzero  eigenvalues  of  and  coincide. 

We  shall  now  briefly  describe  the  relations  between  various  discrete  and 
continuous  Green's  functions.  I/3t  the  discrete  analogue  of  the  logarithmic 
Green's  iunction  be  denoted  momentarily  by  G.  This  discrete  Green's  function  G has 
been  studied  in  groat  detail  in  [d],  [24],  Chapter  3 of  [ 27  ] and  section  4 of  [ 28  ] . It  is 
translational  invariant  so  that  we  may  assume  that  the  second  parameter  is  fixed  at  the 
origin  and  define  G as  a function  of  one  parameter  by 

G(a,b)  = G(P;0)  whore  P-  (a,b),  Oh  (0,0)  . 

Let  G : (a,b)  * R bo  defined  by 

G (a,b)  = (1/2it)  log(a^  + b^)  . 

Let  G , G bo  defined  by 
x’  y 

G^(a,  b)  =■  G(a  + h,  b)  “ C(a,  b)  ; 

Gy(a,b)  = G(a,b  + h)  - G(a,b)  . 

* « 

Let  G , G be  similarly  defined.  It  is  shown  in  section  4 of  [ 28]  that  for  any  nonnegative 
X y 

integers  r and  s,  r > s,  the  following  holds 

( S.  8)  max{  |g  (sh,  rh)  - G (sh,  rh)  I , !g  (sh,  rh)  - G (sh,  rh  I } < (0.  34)r  ^ . 

XX  y y 

(5.9)  G (sh,rh),  G (sh,  rh),  -G  (sh,  rh)  are  always  positive  . 

X y ^y 
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(5. JO) 


-Gyy(sh,rli  - h)  and  G^^(sJi,rh  - h)  - G^^(sh  - h,rh  - h)  are 
always  nonncgative  for  r>s,  r^O  or  sv'O. 

The  following  property  of  G permits  us  to  extend  (S.8)-(5.10)  to  negative  values  of 
r and  s 


(5-11)  G{rh,  sh)  = G(sh,rh)  = G(-sh,rh)  = G(sh,-rh)  . 

Finally,  G satisfies 

1,  r = 0,  s = 0 : 


h A^G(rh,sh)  = 


0,  otherwise  , 


so  that  by  (5.11),  the  following  holds 


(5.12) 


G (0,0)  = G (0,0)  = 1/4  . 

A y 


Let  G‘  be  the  Green's  function  on  the  rectangle  with  zero  Dirichlet  boundary  conditions.  It 

4r 

is  shown  on  p.  515-318  in  [ 11  ] that  G'  and  G differ  only  by  a smooth  function  11  and 

that  B ^ and  G differ  by  a mesh  function  11^^  which  is  an  0(h)  approximation  to  11. 
Using  the  same  technique  of  proof  used  in  [11],  it  is  easy  to  see  that  if  P,  Q are  both 
bounded  away  from  aS  uniformly  in  h,  then  llj^(p,  Q)  is  an  O(h^)  approximation  to 
H(P;Q).  In  what  follows,  wo  shall  denote  B ^ and  G'  by  G and  G unless  stated 
otherwise. 

* 

Lemma  5.3.  aG  /ai'p  is  uniformly  continuous  with  respect  to  both  parameters  P and  Q 

* 

of  G if  both  P and  Q lie  on  a closed  curve  with  continuously  turning  tangent  and  with 
continuous  and  bounded  curvature. 

★ 

Proof.  This  result  is  well  known  if  G is  the  logarithmic  potential.  See  o.g.  [25].  Since 

* 

G'  and  G differs  only  by  a smooth  function,  the  lemma  clearly  follow's. 

6 * 

Lemma  5.4.  Let  P and  Q bo  two  points  in  an,  with  d(P,  Q)  = h , p < 1/2  and  let  P 

h 

* 

and  Q be  their  corresponding  points  on  an.  Let  Op  and  be  the  angles  that  the 

normals  through  P and  Q resriectively  make  with  the  closest  coordinate  axis 
(5.13)  Mp,q)  = 2[aG*/a‘'  J(P*:Q*)  h sec  f o(h^'^P)  . 

p V 
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Kroof . We  shall  only  treat  the  case  for  scheme  l.N.b  when  P has  only  one  neighbor  in 


♦ 


(Cn),  . The  proof  for  all  other  cases  is  almost  identical  and  will  not  be  given.  Let 
h 

W t (C'.f)  . Then 
h 

Q)  = (sec  o^sec  ap)[  2G(W,Q)  - 2(1  - tan  «p)G(P,  Q)  - 2 tan  «pG(S,  Q)]  . 

By  (5.8)  and  the  discussion  following  (5.12), 

G(W,Q)  - G(P,Q)  = G*(W;Q)  - G*(P;Q)  + 0(h^  , 

G(S,Q)  - G(P,Q)  = G*(S;Q)  - G*(P:Q)  ^ O(h^)  . 

Here  6 = min{2,  3(1  - p)}.  Since  the  modulus  of  the  second  partial  derivatives  of  log  d(P,  Q) 

-2 

is  not  greater  than  [d(P,  Q)]  , it  is  easily  veri fied  that  ( 5. 13)  holds. 

Lemma  5.  s.  The  family  of  operators  is  collectively  compact  on  C[0,1], 

Proof.  We  first  show  that  is  collectively  comijact  on  C[0,1].  We  construct  a family 

•k 

of  operators  K on  [ C(  0,  1 ] ] by  the  same  procedure  described  above  with  k(t,  t ) in 
m j 

★ 

(5.1)  replaced  by  k (t,  b)  defined  as  follows 

(5.9)  k*(t,t.)  = yh  sec  «p  [ aG*/a>'](P*_j:P*) 

+ (1  - y)h  sec  <vp  I aGV8vl(P*;P*)  , 

where  the  norma)  derivative  is  taken  with  respect  to  the  first  var  'Ole  and  y is  defined  by 

(5. 2) -(5.  *l).  Let  Ilk  II  = max|k(t,t.)|  and  Ilk  II  = maxllk  ||  . By  lemma  5.3, 

j ^ j ) 

II  K fll  < constant  Ilk  II  Ilf  II  , 

m 

l(K  f)(f)  - (K  0(t')l  < constant  Ilk  - k II  ||f||  . 
m m = t t' 

, * 
i Hence,  is  collectively  compact  on  C[0,1]  by  the  Ascoli-Arzela  theorem.  By 

^ lemma  6.  3, 

„ in  ^ 

II K - K II  < max  ^ |k(t,t  ) - k (t,t  )|  = 0(1)  as  m - 00  . 
fnm  1 1 

-14- 


Therefore  {K  } and  similarly  {K*  } arc  collectively  compact  on  C[0,ll.  The  theorem 
m m 

easily  follows. 


Lemma  fi.  6.  K_f  — K Kf  for  each  ft  C[0, 1)  where  K is  the  compact  integral  operator 


defined  by 
(5.10) 


(Kf)(t  ) = 2 / ( 8G  AvjJf  ds  , 

an 


where  P = (<!)(tp),  4'(tp)). 

Proof.  Let  {Q.  s Pj^^,  j = l,...,n}  be  a subset  of  352^  that  is  chosen  as  follows.  The 

are  strictly  increasing  as  j ranges  from  1 to  m and  between 

, * 

and  2vh.  Let  0^  be  the  corresponding  points  on  all.  It  is  easily  seen  that  as  h — 0, 

(5.11)  (Kf)(tp)  =2^1  aG*/av'](P:0*)d(Q*,Q*^j)f(t.^.j)  + 0(1)  . 

* * 

The  ^(^j  I on  the  right  hand  side  of  (5. 11)  can  be  replaced  by  [ i(j  + 1)  - i(j)  ] h sec  Cq 

without  affecting  the  0(1)  nature  of  the  remaining  term.  It  then  easily  follows  from  lemma  5.4 

that  K f Kf  and  similarly  K'  f -►  for  each  f c C[0,1). 

m m 

Proof  of  Theorem  5.1.  By  lemmas  5.1,  5.  5 and  5.6,  we  see  that  Theorem  5.1  holds  for  the 

singular  values  of  K . The  theorem  then  follows  because  of  lemma  5.2. 
m 


Theorem  5.2. 


0. 251  < B,  B,  < 7 . 291  for  scheme  I.  N.  a ; 
— h h ~ 

0. 2 51  < bJb  < HI  for  scheme  I.  N . b . 

~ h h ~ 


Proof.  Below  we  give  the  proof  for  scheme  l.N.a.  Details  of  the  proof  for  scheme  I.N.b  may 
be  found  in  section  5 of  [28],  where  it  is  shown  that  the  B^^  matrix  for  sclicmes  l.N.a  and 
I.N.b  are  essentially  the  same  as  that  of  schemes  I.  a and  I.b  considered  in  [28]  respec- 
tively. We  shall  first  prove  that  the  follov/ing  holds  for  scheme  l.N.a 


(5.14) 

The  following  lemma  is  well  known. 
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Lemma  5.7.  Let  the  symmetric  part  of  a matrix  A satisfy 

(A  + A^)/Z  >61,  6 > 0 . 

Then 

T 2 
A A > 6 I . 

T 

Let  B denote  + B.  . We  shall  show  that 
s h h 

(5.16)  min  {B  (P,  P)  - ^ |b  (P,0)I)>1 

PcaQj^  QcaQj^,Q^p 

so  that  (5.14)  holds  because  of  lemma  5.7  and  a well  known  Gerschgorin  theorem. 

Let  P € 8Q  . Assume  that  the  local  orientation  of  the  boundary  near  P is  such  that 
h 

for  any  point  P'  € in  a neighborhood  of  P,  either  W and  N',  the  western  and 

northern  neighbors  of  P',  are  both  in  (Ct2)j^  or  W'  alone  is  in  (CJ7)^.  Let  E and  S 
denote  the  eastern  and  southern  neighbors  of  P respectively.  Then 
(5.16)  B^(P,  Q)  = (sec  Og/sec  ap)(G^(W:Q)  - (tan  a^)  G^(S;Q) 

+ G (N;Q)  - (tan  a ) G (E;Q)]  . 

X * y 

Here  the  subscripts  x and  y denote  the  forward  differences  in  tlie  x^  and  x^  directions 
respectively  taken  with  respect  to  the  first  parameter  of  G.  Because  of  the  band  structure 
of  B^,  we  may  assume  without  loss  of  generality  that  G is  the  discrete  analogue  of  the 
logarithmic  Green's  function.  The  error  resulting  from  this  assumption  is  less  than  a constant 
times 

Assume  that  tan  is  bounded  away  from  0 and  1.  Then  for  h sufficiently  small, 

(an  ) , a subset  of  90^^  which  contains  a Vh  neighborhood  of  P can  be  partitioned 

loc 

into  blocks  as  follows.  Let 

Ijj  = {(0,h),  ...,(0,Mjh))  , 

= {(kh,  M^h  + h) (kh,  Mj.^jh)},  k = 1,  . . . , Kj  , 

I_l^  = {(-kh,  • • • .("kh,  -M_^h  + h)},  k = 1, . . .,  . 
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Then 


(an  ) = U 1 : p » I . 

n , 1 k 0 

loc  *^""^2 

Let  and  ni_|^  be  the  number  of  points  in  the  nonempty  sets  ^ and  I ^ respec- 
tively. Then  M = m,  + • • • + m,  : M , = m , + •••+  m , : = 0.  Let  P denote  the 

k 1 k -k  -1  -k  0 ) 

point  with  y-coordinatc  jh.  By  (fi.8),  (“j-lB)  and  the  smoothness  of  an,  it  is  easily 

verified  that  I®  (^iQ)!  will  remain  essentially  unchanged  if  tan  is  replaced 

Q^P  ® ^ 

throughout  by  tan  np.  Let  a ^ tan  dp.  Let  G(i,  j)  b G(ih,jh).  Let  P e P, . Then 

(5.17)  B (P,  P)  = 2 + 2(1  - a)G  (0,0)  ; 

s xy 

B^(l’.P^)  = (1  - a)[Gy^(l,  |j  - i|  - 1)  - Gyy(0,  Ij  - i|  - 1)],  Pj  ‘ Iq  . 

Hence,  if  P.  c ]^,  then 

(5.18)  - a)lG^(0,0)  - G^(0,  - i)  - G^(0,1)  + G^(0,  Mj  1 1 - i)] 

(5.19)  ^ B (P,  P ) = (1  - a)|G  (0,0)  - G (0,i  - 1)  - G (0,1)  + G^(0,i)]  . 

. _ . ® J X X X X 

1>) 

Similarly,  if  P^  t I^,  k = l,2,...,Kj,  then 

Bg(P,P.)  = (1  - a)[G^y(k,j  - i)  4 G^y(k  - l,j  - i - 1)1  , 


4G^(k-l,M^-l)-G^(k-l,M^^j-l)J  . 


On  the  other  hand,  if  P t I , k = 1,2,...,K  then 

“J  ” N C 


Bg(P.  P_j)  = (1  - a)[  G^  Jk,  i 4 j)  4 G^  Jk  - 1,  i 4 j - 1)  ] , 


Z B,(f.  I'.,)  ■ -0  - .)l  G^(k,  * .)  - G^(k,  M_^  t 1) 


-J  -k 


4 G (k  - 1,  M . 4 I - 1)  - G (k  - 1,  M 4 i - 1)  J 


■J 


It  is  easily  verified  from  (5.9)  and  (5.10)  that 

Bg(P,Pj)<0,  P.  c k = ±l,±2,-.- 
(5.20)  B^(P,P.)>0,  P ‘ • 

Since  we  only  want  to  obtain  an  upper  bound  for  V |b  (P,  P.)  I,  we  may  assume  without 

loss  of  generality  that  B (P,  P.)  * 0 iff  P^  t (an,  ) . Then,  by  summing  P,  over  j > M 

® ^ ^ loc  J ^ 

- Yj  Pj)  = (1  - a)[  + 1 - i)  + G^(0,  - i) 


+ y {G  (k,  M,  + 1 - i)  + G (k  - l.M.  - i)} 

^ XX  ’ k+1  ifv'  * k4! 


By  (5.9) 


>Mg  (k  + 1;M,  - i)  + G (k,  M,  - 1 - i)l  < G (2,  - i)  + G (1,  M,,  - i - 1)  . 

t-j  '■  vv'  k+1  yy'  ’ k+1  V ’ 2 V ’ 2 


y'  ’ 2 


y'  ’ 2 


Hence,  the  following  holds  when  P^  is  summed  over  all  j > 

(5.21)  ^ Ib^(P,P  )l  < (1  - a)[G^(l,  Mj  + 1 - i)  ^ G^(0,  - i) 


+ Gy(2,M^  - i)  H Gy(l,M^  - i - D)  . 

Similarly,  the  following  holds  when  P^  is  summed  over  all  j < 0 

(5.22)  Y < 0-a)(G^(l,i)  +G^(0,i-1)  + G^(2,  M_j  + i - 1)  + G^d,  + i - 2) ] 

By  (5.17)-(5.22), 

(5.23)  B (P,P)-  Y Ib  (P,P  1|  < 2 - 4(1  - a)[G  (0,0)  - C (0,1))  - (1  - a)  H(i)  , 

“ Pyp  s J X X 


H(i)  = G^(0,  i)  + G^(l,  1)  + G^(0,  + 1 - 1)  + G^(l,  Mj  + 1 - 1) 


+ Gy(l,  M_j  + 1 - 2)  + Gy(2,  M_j  + 1 - 1)  + G (1,  - I - 1)  + 0^(2,  M^,  - i)  . 
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It  is  easily  seen  from  the  table  on  p.  292  of  [21]  that  I!(i)  attains  its  maximum  at 


I 

( 


either  1 = 1 or  1 = M^.  It  is  then  easily  verified  with  the  aid  of  the  above  mentioned 
table  that 

(5.24)  4aG  (0,0)  + 4(1  - a)G  (0,1)  > (1  - a)H(i),  i = l,...,M  . 

XX  1 

By  (5.  23)  and  (5.  24),  we  see  that  (5. 15)  holds  for  our  choice  of  P. 

The  proof  for  other  choices  of  P is  almost  identical  and  will  not  bo  completed.  Thus 
we  complete  the  proof  for  (5.  6).  We  now  proceed  to  prove 

(5.25)  < (7.29)1  . 

n n 

We  shall  show  that  the  following  holds  * 

(5.26)  max(2^  I K.  (P,  Q)  I , Ib![(P,  Q)  I } < 2.  7 . 

g " Q 

We  assume  that  the  local  configuration  of  points  in  9f2,  in  a s/h  neighborhood  of  P is 

h 

★ 

Similar  to  that  in  idU.)  described  earlier  in  this  section.  Let  B,  (P,  Q)  be  defined  by 
h , h ’ 

* j * 

equation  (5.8)  with  G replaced  by  its  continuous  analogue  G . Let  B^^  (P,  O)  be 
similarly  defined.  We  shall  first  show  that 

(5.27)  Y Ib*(P,Q)I  <0.7  . 

lygT>3 

Lot  r,  r*,  r^,  r',  r„  and  r^  denote  d(P,  0),  d(w,  Q),  d(S,  Q),  d(NE,  Q),  d(N,Q)  and 
d(E,  Q)  resfiectively.  We  again  assume  that  tan  Oj,  s tan  s a for  all  Q such  that 
B^(P,  C?)/0.  Then 

2ir  • Bj^(P,  Q)  = -log(r  /r)  + a log(rj/r) 

-log(r,/r')  + a log(r2/r')  . 


Let  the  coordinates  of  P,  Q be 

*2  2 ^.22 
r - r = 2hxQ  + h ; r^  - 


(0,0)  and  (^Qi  Vq)  respectively.  We  have 

•2  2 2 2 2 2 '2 
r = 2hXp  - h ; fj  - r = 2hyp  + h ; r^  - r = 
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Hence, 


4"B^(P,  Q)  - (-2x^h  + 2aygh)(l/r^  ^ 1/r  ^)  + (a  - l)h^(l/r^  - 1/r  ^)  + Rg  , 

where 

IRqI  < t(2Xg  + h)^  + a(2yg  + h)^)hV2r‘’  4 l(2Xg  - h)^  4 o(2yg  - h)^]hV2r  . 

It  is  easily  verified  that  if  P,  Q are  two  jx)ints  on  an  with  d(P,  O)  < h^,  0 < y < 1, 
and  ta  is  the  tangent  at  0 to  3J2,  then 

d(P,  to)  < (k  4 0(l))h^^  , 

where  k is  the  maximum  absolute  value  of  the  curvatures  of  an.  Hence,  for  all  practical 
purposes  we  may  assume  k = 0.  Then 

Xg  = ayg  4 eh,  |e  1 < 1 . 

Hence, 

^ hrr  . B*(P,Q)I  < P(r,r-)  , 

lygl>2 

where 

F(r,r')^  21  {2h^{l/r^  4 1/r'^)  4 (1  - a)|l/r^  - 1/r’^l 

lygl>2 

22  224  M 444  *4 

4 (2Xgh  4 2aygh  )(l/r  4 1/r  ) 4 (1  4 a)(h  /2r  4 h /2r  ) 

4 2(Xg  4 ayg)(h  V2r‘*  - hV2r  "*)}  . 

Talcing  into  account  the  local  configuration  of  points  in  a Vh  neighborhood  of  P for 

a certain  a b tan  o^,  we  easily  see  that  for  all  0 < a < 1, 

1 ^ - 1)^  4 kV(lc  4 1)"*}  . 

|yg|>2  ^ k=3 

Hence, 
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r(r,r')<  2]  (8A^  + (3  - a)/(k  - 1)^  + (1  + a)/(k  + 1)^ 

k=  3 

+ 2k^(l/(k  - 1)‘‘  + l/(k  + D^*) 

+ lA**  + l/2(k  - o'*  + l/2(k  + I)'’ 

+ 2*s/2  (lA^  + l/(k  + 1)^} 


This  completes  the  proof  of  inequality  (5.27).  By  (5.17), 

(^•28)  IBjj(P,P)|  <1  . 

It  is  easily  verified  that  £ Ib^(P,  Q)  1,  P ^ Q,  |yQ|<2  attains  its  maximum  when  a=0 
and  P has  two  neighbors  in  (^$2)^^.  A simple  calculation  with  the  aid  of  table  I in  [ 4j 
or  table  II  in  [24]  shows  that  when  summed  over  all  Q with  0 < lypl  < 2 

(5-29)  2|b^(P,Q)|  <0.96  . 

By  theorem  4.3  of  ( 28)  and  table  I of  ( 4] , we  see  that 

(5-30)  X,  Ib.(P,Q)  - B*(P,Q)|  <0.04  . 

Iy„l>3  h h - 


By  (5.27)-(5.  30), 

X Ib.  (P,Q)I  < 2.7  . 

Qcan.  ^ 
n 

Since  the  same  inequality  holds  by  a similar  argument  when  B in  (5.  31)  is  replaced  by 
T 

Bj^,  we  easily  see  that  (5.  26)  and  hence  (5.  25)  holds. 

.Umma.3.8.  If  D = A + B,  where  A and  B are  arbitrary  matrices  with  singular  values 

"l  - “’2  - ” ’ - “n  - ° Pj  ^ - " ■ - ® respectively  and  6>6  >-*->6  >0 

are  the  singular  values  of  D,  then 

*1+J+1  - "i+1  ^j+l»  ^ positive  integers  . 

Proof.  See  e.g.  ex.  28  on  p.  89  of  [20]. 
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•Theorem  5.  3.  Let  dj  and  d^  be  the  spectral  bounds  of  B^.  Then  given  e > 0,  there 
exists  a positive  integer  N independent  of  h such  that  all  except  N singular  values 
of  C lie  in  [dj  - e,d^  + c). 

Proof.  Let  C.  K.  and  B.  be  matrices  A,  B and  D respectivelv  in  lemma  5.8.  Bv 
h h 

theorem  5.1,  given  e > 0,  there  exists  a positive  integer  p such  that  for  all  h > 0, 

Vj+1^^-  ^ 

Since 

d,  < 6 < a .-ip  j = 0,1,2,  . ..  , 

1“  m m-p'j  p-fj+r  > > « > 

only  the  last  p -f  1 singular  values  of  C may  lie  to  the  left  of  d^  - e.  Similarly,  by 

letting  B.  , -K,  and  C be  the  matrices  A,  B and  D respectively  in  lemma  5.8,  we 
n n 

see  that  only  the  first  p + 1 singular  values  of  C can  lie  to  the  right  of  d^  + e. 

Theorem  5.  4.  Let  I|m1|  and  R(M)  denote  the  spectral  norm  and  range  of  a matrix  M 

respectively.  Let  A,,  be  the  same  as  in  (3.2).  Let  A = aJ.A,,  |R(Ay.)  and 
11  s 11  11  11 

C = c'^cIr(C^)  be  the  restrictions  of  A^A  and  c'^^C  to  R(aJ)  and  R(C^)  rospec- 
S ii  11  11 


tively.  Then 


Ic-^ll  < (Ia-MI  IIbII^IId-^II^ 

s s 


where  B and  D are  the  same  as  in  equation  (3.  5). 

T T 

Proof.  Let  v = UU  v t R(A).  It  is  shown  in  section  3 that  U v < R(C)  so  that  we  may 

make  the  Ansatz  u = B ^UDp  for  the  solution  of  Au  = v and  solve  the  alternative  form 

T 

of  capacitance  matrix  equation  Cp  = U v for  p.  Ixjt  Pj  be  the  eigenvector  corresponding 

T 

to  the  smallest  eigenvalue  of  the  positive  definite  matrix  C^.  Let  Cpj  = U v^,  where 


Vj  = UU  Vj.  Then 
(5.32) 


Ip.ll''^  llc"MllluV||" 


Here  l|v||  denotes  the  Euclidean  norm  of  a vector  v.  Let  = B UDpj.  Let  U and 
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U*  be  the  extension  operators  from  Sf.  U an.  to  all  mesh  points  and  from  n.  to 

h n n 

* ♦T  * * 

n.  U an^  respectively  that  are  defined  the  same  way  as  U.  Let  u,  = U u , = u_  + u, 
h h 1 1 3 

* ♦ T . 

where  t the  null  space  of  ‘ R(A^j).  Because  of  the  reducible  structure 

*T 

of  A,  U VC  N(Ajj)  if  V < N(A).  Hence  we  may  write  *^3>  where  u^  c N(A), 

*j  * *-p  * 

U *^2  ” *^2  ^ *^3  ~ *^3'  ip  * 0 >■  N(C).  It  is  readily  verified  that 

B ^UDy*  ^ 0 c N(A).  Since  both  N(A)  and  N(C)  are  of  dimension  one,  we  can  choose  v'j 
such  that  B = u^.  Then  u^  = U ^B  V U'D(p^  - <p^).  Clearly,  ||  - <,"^11  > lipjll- 

Hence, 


(5.33) 


* *J 

On  the  other  hand,  ^12*^3  ~ ^ ''j-  Hence, 

(5.34)  ||u!||2  < llA-'llllv  II'  . 

3 S 1 


The  theorem  then  easily  follows  from  (5.  32)-(5.  34). 
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§ 6.  Rate  of  convergence  of  tlio  coniuciiitp  erradipnt  iterations 

Let  b denote  the  right  hand  side  of  the  capacitance  matrix  equation.  We  arc  concemoc 
T T 

with  solving  C Cp  = C b by  the  conjugate  gradient  method.  It  is  shown  in  [21]  that  the 
conjugate  gradient  method  gives  the  solution 

p = C'^b  + (I  - P)pj^  . 

Here  is  the  generalized  inverse  of  C.  For  any  m-vector  b,  C^b  is  the  unique  least 

square  solution  of  Cp  = b that  is  of  the  minimum  Euclidean  norm.  P is  the  orhtogonal 

T 

projection  of  onto  R(C  ) and  p^  is  the  initial  guess.  We  shall  assume  thioughout 

T + T -1 

this  section  that  p^  e R(C  ) so  that  p = C b.  Let  v = C b.  Then  p = v.  More- 
over, from  (6. 1)  and  (6.  2)  below  we  see  that  all  the  relevant  vectors  generated  in  the 

T 

conjugate  gradient  process  are  in  R(C  ).  Hence,  the  original  problem  is  reduced  to  solving 

Qp  = V by  the  conjugate  gradient  method  where  Q ^ C is  a positive  definite  symmetric 

s 

matrix. 

We  now  briefly  describe  the  conjugate  gradient  method.  Seee.g.  [9],  [M],  [19], 

[16]  and  [ 23]  for  details.  Let  p^^  = -g^^  = v - Qpq.  The  conjugate  gradient  process 
generates  a sequence  of  vectors  Pj^  approximating  the  solution  p by 

Pk+i" 

Pk+1  = -^k+l^t^k+l^f'k/pk'^pk^pk- 

T 

whore  gj^  = Qx^  - v.  The  Pj^  are  0-conjugate,  i.c.  p^  Op.  = 0,  i ^ j.  The  p^ 

T T 

minimizes  the  quadratic  form  (l/2)w  Ow  - v w on  the  linear  variety  p^  + where 
Is  the  subspace  spanned  by  {p^,  Pj,  . . . , P|^  j}.  The  iterates  pj^  satisfy 

(6.3)  Pk  = PQ  ^ 

where  j Is  a polynomial  of  degree  k - 1.  It  Is  shown  In  [ 23]  that  among  all  iterative 
methods  that  satisfy  (6.  3),  the  conjugate  gradient  method  is  optimal  in  the  sense  that 

'<1 
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I, 

4 

i 

I ■ 


E(P|,)  = (1/2)(P^  - p)’'q(p^.  - p) 

is  minimal.  It  then  easily  follows  that  if  \.  are  the  eigenvalues  of  Q,  then 

for  any  choice  of  a polynomial  of  degree  k - 1.  Lot  Z : (a,b)  ^ R,  where  (a,  b) 

are  ordered  pairs  of  positive  numbers  a and  b be  defined  by 

Z(a,b)  = [(I  - N^a)/(1  + Va)l^  . 

It  is  known  (sec  e. g . [ 9)  or  ( 23])  that  we  can  select  ^(k)  such  that 

max|l+\.P  (k  .)  I < 2Z(k,  k)  , 

» 1 1 1 

’^i 

where  k is  the  spectral  condition  number  of  Q.  On  the  other  hand,  suppose  all  except 
N eigenvalues  of  Q lie  in  the  interval  (c^.c^j.  Lot  k . , i = l,  . . . , N be  the  excep- 
tional eigenvalues,  let  Then  as  before  we  can  choose  ) •‘5uch  that 

< 2((i  - k^K,)/(l  + 

c,  <k.  <c,  ' ^ ^ 

1 i~~  2 

Choose  such  that 

Then, 

N 

(6.6)  max  |l  +k^pj^_^(k^)  I < 2Z(Kj,k  - N)  max  { TT  (l-k/k  |}  . 

Cj<k<c2  i = l ‘ 

By  (6.4)-(6.6)  and  theorem  S.  3,  we  easily  have  the  following  theorem. 

Theorem  6.1.  Let  k and  k,  be  the  spectral  condition  numbers  of  C and  resoec 

^ s h h 

lively.  Then  given  c > 0,  there  exists  a positive  integer  N independent  of  k and  h 
such  that 
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E(p^)/E(Po)lmin{4Z(K,2k),-lZ(Ki  - e,  2k  - 2N)x(\)}  . 


Horo  x( \)  - max  7T  whore  i = arc  the  N eigenvalues 

Cj  < k <c^  i = l * '■ 

of  C that  lie  outside  [ c.  - e,c  + e].  Here  c,  and  c,  arc  the  spectral  bounds 

S A ^ 1 ^ 

of  . 

h h 

Corollary  6.1.  The  number  of  iterations  needed  to  reduce  E(P|^)/E(pq)  to  a given  accuracy’ 
can  grow  no  faster  than  constant  log  m as  h — 0. 

Proof.  By  theorem  5.-i,  the  smallest  eigenvalue  of  is  larger  than  constant  - m The 

corollary  is  therefore  an  immediate  consequence  of  theorems  5.2,  5.  3 and  6.1. 


. § 7 . Survey  of  previous  work  on  capacitance  matrix  methods 

We  give  here  only  a brief  survey  of  previous  work  on  capacitance  matrix  methods. 

See  also  section  8 of  [26]  for  more  references  and  more  details  on  some  of  the  references 
mentioned  here.  C.  W.  Hockney  gave  a brief  description  of  a method  of  this  type  in  [18]. 
He  credited  Oscar  Buncnian  for  the  idea. 

The  papers  [7  ] and  [13]  by  Buzbee,  Dorr,  George,  and  Golub  and  George  respectively 
used  the  same  Ansatz 

u = B”*v  -1  B~^UDp 

as  is  used  in  our  algorithm  to  treat  tJio  Dirichlet  problem.  It  is  then  shown  exporirncntally 
in  [26]  and  theoretically  in  [ 27]  that  the  resulting  capacitance  matrices  C are  ill- 
conditioned  and  that  tlie  singular  values  of  C cluster  around  zero.  The  conjugate  gradient 
method  was  used  in  [13]  to  solve  the  capacitance  matrix  equations  using  an  iterative 
imbedding  technique  similar  to  that  mentioned  in  section  -1  of  this  work.  The  number  of 
iterations  used  to  achieve  a given  accuracy  are  proportional  to  the  square  root  of  m,  the 
order  of  C.  The  regions  considered  in  [7]  a.''e  of  a rather  simple  type.  The  matrices  C 
are  jxisitive  definite  syiinietric  and  the  Cholesky  method  is  used  to  factorize  C.  The 
j numerical  result.s  aie  obtained  on  a C13C6600  and  a gain  in  speed  of  a factor  three  is 

j reiK)rtod  in  [ S]  lor  runs  on  GDC  7 600. 

The  paper  [ 26]  by  I’roskurowski  and  Widlund  is  probably  the  first  one  that  exploits 
I the  similarity  between  the  classical  potential  theory  and  the  capacitance  matrix  method. 

I It  is  shown  experimentally  tlicre  that  by  making  the  correct  Ansatz  guided  by  the  classical 

' potential  theoiy  the  capacitance  matrix  method  becomes  a well  pxjsed  problem.  The  matrices 

, C for  many  test  regions  arc  uniformly  well  conditions  in  the  spectral  norm  and  the 

convergence  of  the  conjugate  gradient  iterations  for  tl.  se  regions  appears  to  be  Independent 
of  the  mesh  sizes.  It  Is  then  shown  theoretically  in  { 26]  that  (or  a large  class  of  domains 

I 
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and  some  special  schemes  of  approximating  boundary  conditions,  the  above  observation 
is  always  valid.  This  work  is  mainly  an  outgrowth  of  [ 26].  It  extends  some  of  the 
theoretical  estimates  in  [26]  to  all  bounded  domains  with  smooth  boundaries. 


I 


§ 8 . Numerical  cxucrimonts 

To  Illustrate  the  effectiveness  of  our  algorithm,  wo  have  used  linear  polynomials  as  test 
functions.  Truncation  errors  are  not  present  and  the  right  hand  side  of  equation  (8.1)  is  always 
consistent 

(8.1)  A^jU  = U*’’v  . 

it 

Let  u and  u be  respectively  the  exact  and  numerical  solutions  of  (8.1).  Let 

E £ II 6 (u  - u )|l  + life  (u  - u )ll  , whore  fe  and  fe  denote  the  undivided  forward 

max  Xj  x^  “ Xj  x^ 

differences  in  the  x^  and  x^  directions  respectively.  The  domain  Q is  an  ellipse  with  the  ratio 
of  half  axes  equal  to  y and  the  test  function  u satisfies  u(x)  r x^.  The  following  is  a table  of 
numerical  results  obtained  by  test  runs  on  the  Univac  1110  at  MACC,  University  of  Wisconsin, 
Madison. 


Table  1 


No.  of  iterations 

V 

m 

Norm  of  C.G.  Residual 

E (approx.) 

max  * * 

4 

1 

36 

. fel720fe6-03 

-04 

4 

1 

76 

. fe40990fe-03 

-04 

4 

1 

108 

.7231133-03 

-0  1 

4 

0.7 

32 

. 3189fel0-02 

. 2-03 

4 

0.7 

64 

. 86fe8407-03 

.2-03 

4 

0.7 

92 

. l2667fe3-02 

.2-03 

4 

0.  fe 

60 

. fe? 68389-02 

. 2-03 

4 

0.  fe 

84 

. 3497684-02 

. 2-03 

7 

1 

108 

. 1820106-04 

-06 

7 

0.7 

92 

. 1 37201fe-04 

-06 

7 

0.  fe 

84 

. 3218270-01 

-06 

Scheme  l.N.a  is  used  to  obtain  results  listed  in  Table  1.  Typically  it  will  take  one  or  two 
more  iterations  to  achieve  similar  accuracies  if  scheme  I.N.b  is  used.  The  norm  of  C.G.  residual 
given  in  the  fourth  column  of  Table  I is  the  L^  norm  of  the  conjugate  gradient  residuals  divided 
by  the  square  root  of  the  number  of  mesh  points  Inside  n. 
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• § 9 . Conclusions 

Since  it  takes  two  fast  Poisson  solvers  to  complete  each  conjugate  gradient  iteration, 

the  total  operation  counts  of  the  algorithm  are  approximately  ten  or  eleven  times  that  of 

a fast  Poisson  solver  for  the  Laplace  or  Poisson  equation  respectively.  It  is  reported 

in  (2)  that  the  operation  c^ounts  of  a fast  Poisson  solver  can  be  reduced  to  O(n^)  if 

the  fast  Pourier  transform  methods  are  combined  with  k cyclic  reduction  methods  if  k 

is  projxjrtional  to  log^n.  It  is,  however,  more  realistic  to  say  liiat  the  operation  counts 

2 

of  our  algorithm  are  proportionally  to  n log^n  in  the  experiments  carried  out  so  far. 

2 2 

Our  theoretical  estimate  of  a constant  times  n (log^n)  is  perhaps  too  conservative.  It 

2 

is  shown  in  [ 27]  and  [ 28)  that  theoretical  estimates  of  constant  • n (loy^n)  can  bo 
obtaitied  for  a special  clas.s  of  do.mains  in  some  cases  although  it  is  still  an  open  question 
whether  this  is  true  in  general. 
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of  the  algorithm  docs  not  exceed  constant  times  n-(log  n).  (n  = l/h)  for  any 
bounded  domain  with  sufficiently  smooth  boundary. 
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